

%% Maximize the Likelihood function
%% Changes for the nopq version -> N=1 (no price states), so total states = X*X
%%				-> p,q do not enter into U
%%				-> impchoicestate is just prev x, future x.

clear

cd ../china/data
load hs10_indlist.raw;
indN=length(hs10_indlist);

for w=1:50


folder=sprintf('../china/estimation/hs10/regular/raw_%d',hs10_indlist(w,1))


Name='Ryan is Great';

tic
addpath('/apps/tomlab');
startup

cd (folder)

load data2_q_regular_nopq.mat
N=1;
nExp = X;
nImp = M;
%S=N*nExp*nExp;
S=nExp*nExp;
delta=0.975;

failnum=0;


% tom objects include beta_p, beta_x and EV
%beta_p=tom('beta_p',1,1);
beta_x=tom('beta_x',1,1);
beta_c=tom('beta_c',1,1);
%xi_tilde=tom('xi_tilde',1,1);
EV=tom('EV',S,1);


U=beta_x*switched + beta_c*city_switched ;

% EV is a tomSym object, so cannot write:
% for i=1:S CSVF(i,1) = U(i,1) + delta*EV(i,1)

CSVF= U + delta*EV;

index=ones(nExp,1)';
for i=1:nExp
    index(1,i)=N*nExp*(i-1)+1;
end

index=repmat(index,N,1);
index=index(:);
index=repmat(index,nExp,1);

% index is a vector that tells us what the first number of EV calculation
% is for each element of EV.

%TransProbs=rand(N*N*nExp*nExp); Too big!

A=zeros(nExp,N,S);

for i=1:S
    B = (index(i,1):index(i,1)+N*nExp-1)';
    B = reshape(B,N,nExp);
    A(:,:,i) = B';
end
% Each A is a (nExp x N) matrix arranged so as to calculate EV
% For example EV(1) is a double sum: sum of  sum of CSVF(1),(6),(11);
% sum of CSVF(2),(7),(12); etc.

Q2 = zeros(nExp, S);
for j=1:nExp
    Q = A(:,:,N*(j-1)+1);
    Q2 (:, (j-1)*nExp*N +1: j*nExp*N ) = repmat(Q,1,nExp);
end
% Q2 is the matrix for adding up to calculate the denominator of P



 P = exp(CSVF) ./ sum( exp( CSVF(Q2) ) )' ;



% New Constraint Definition:


W=tomArray(CSVF,[X X]);
Y = log( sum( exp( W),1 ));
Y = repmat(Y,X,1);
Y = Y(:);

%W=repmat(W,N,1);
%W=W(:);
%W=repmat(W,X,1);
%W=tomArray(W,[N X N*X*X]);
%Y=  log( sum( exp( W),2) );
%TransProbs=reshape(TransProbs,1,X*X);
%Z=sum( Y.*TransProbs,1);
%Z=Z(:);




objective = sum ( log ( P ( impchoicestate) ) );

constraints = {EV==Y};


%load beta.mat
%load EV.mat
%EV_guess=mean(EV_sol);
%EV_guess=repmat(EV_guess,S,1);
EV_guess=ones(S,1);

guess = struct('beta_x',-2,'beta_c',-2,'EV',EV_guess);
options = struct;
options.name = Name;
options.prilev=3;
options.use_H=false;
options.use_d2c=false;

Prob = sym2prob('minlp',-objective,constraints,guess,options);

%Prob.DUNDEE.optPar(20) = 1;
%Prob.optParam.IterPrint = 1;
Prob.KNITRO.options.Feastol_Abs=0.0001;
Prob.KNITRO.options.Opttol_Abs=0.0001;
Prob.PriLevOpt=3;

Result = tomRun('knitro',Prob,2);
toc


b=Result.x_k;
beta_c_sol=b(S+1);
beta_x_sol=b(S+2);
EV_sol=b(1:S);

[beta_x_sol beta_c_sol]

save ('beta_regular_nopq', 'beta_x_sol','beta_c_sol'); 
save ('EV_regular_nopq','EV_sol');

w

cd ../china/estimation/hs10

end

diary off